import string
import copy
import scipy
import Tkinter, tkFileDialog
import numpy as np
import pandas as pd
from scipy.spatial.distance import cdist
import matplotlib.pyplot as plt
import os
import sys
import datetime
import matplotlib.units
import re
import glob
from numba import jit,int32
import time
sys.path.append(os.path.abspath("C:\Users\Scherer Lab E\Documents\GitHub\Python_Data_Analysis"))
import common_functions
%matplotlib inline

C:\Anaconda\lib\site-packages\skimage\viewer\utils\core.py:10: UserWarning: Recommended matplotlib backend is `Agg` for full skimage.viewer functionality.
  warn("Recommended matplotlib backend is `Agg` for full "

import itertools

def plot_spiff_different_nn_separations(nn_df, nn_1st, nn_2nd, x_key='x pos', y_key='y pos'):
    nn_combinations = list(itertools.product(nn_1st, nn_2nd))
    len_nn_1st = len(nn_1st)
    len_nn_2nd = len(nn_2nd)
    pivot_table_nn = nn_df.pivot_table(index=['frame','track id'], columns='nn_num', values='nn_dist')
    #print pivot_table_nn.iloc[:,1].mean()
    #plt.hist(pivot_table_nn.iloc[:,1], bins=30)
    #plt.show()
    #plt.plot(pivot_table_nn.loc[:,1], pivot_table_nn.loc[:,2], 'o', alpha=0.05)
    #plt.show()
    plt.figure(figsize=[6*len_nn_2nd, 5*len_nn_1st])
    for num, nn_lim in enumerate(nn_combinations):
        valid_1st_nn = pivot_table_nn.loc[:,1] < nn_lim[0]
        valid_2nd_nn = pivot_table_nn.loc[:,2] > nn_lim[1]
        valid_particles = valid_1st_nn & valid_2nd_nn
        #print nn_lim,valid_particles.sum(), valid_1st_nn.sum(), valid_2nd_nn.sum()
        #print pd.concat([pivot_table_nn, valid_1st_nn, valid_2nd_nn, valid_particles], axis=1)
        temp_nn_df = nn_df.set_index(['frame', 'track id']).ix[valid_particles[valid_particles].index].reset_index()
        
        plt.subplot(len_nn_1st, len_nn_2nd, num+1)
        plt.hist2d(temp_nn_df['x pos']%1, temp_nn_df['y pos']%1, bins=50)
        plt.axis('equal')
        cbar = plt.colorbar()
        cbar.set_label('counts')
        plt.minorticks_on()
        plt.title('1st NN < '+str(nn_lim[0])+' pixels \n 2nd NN > '+str(nn_lim[1])+' pixels')
    plt.tight_layout()

Other Data Set

from IPython.display import Image
os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\SPIFF Manuscript\Data")
Image(filename="06191317_example_image.png")

#os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\Dynamical Phase Transition\Data\SPIFF")
Image(filename="06191317_particle_closeup.png")

Import the Mosaic Data for Mov_06191306 for different windows

Used windows with radius 3, 4, 5, 6, and 7 pixels.

os.chdir("E:\Pat's Projects\Constant Dragging\Exp06191301\Mov_06191306 - 200fps 350x2560")

mosaic_list = []
for traj_data in glob.glob('Mov_061913*'):
    rad_size = int(traj_data[-5])
    df_mos = pd.read_csv(traj_data, sep='\t')
    df_mos = df_mos.rename(columns={'x': 'x pos', 'y': 'y pos', 'Trajectory': 'track id', 'Frame': 'frame'})
    mosaic_list.append([df_mos, rad_size])

Re-Link Particles with Trackpy

import trackpy as tp

'''Link the image data for each window'''
tracked_data = []
for frame, window in mosaic_list:
    tracked = tp.link_df(frame, 10, memory=1, pos_columns=['x pos', 'y pos'])
    tracked_data.append([tracked, window])

Frame 974: 4284 trajectories present

'''Checkpoint after linking'''
os.chdir("C:\Users\Scherer Lab E\Documents\IPython Notebooks\SPIFF Manuscript\Data")
import cPickle
pik = open('Ana_16042401_linked_positions.pkl', 'w')
cPickle.dump(tracked_data, pik)
pik.close()

import cPickle
pik = open('Ana_16042401_linked_positions.pkl', 'r')
tracked_data = cPickle.load(pik)
pik.close()

def disp_calc_en(grp, tau):
    if len(grp) < tau:
        return pd.Series([np.nan])
    #print grp[['x pos', 'y pos']],grp.shift(tau, axis=0)[['x pos', 'y pos']]
    disp = ((grp[['x pos', 'y pos']] - grp.shift(tau, axis=0)[['x pos', 'y pos']])**2).sum(skipna=False,axis=1)
    #if grp.name > 5000:
        #print grp.name
    return disp

def calc_en_msd(df, percent=0.3):
    en_msd = []
    pivot = pd.pivot_table(df, values=['x pos','y pos'], index='frame', columns='track id')
    for tau in range(int(max(df['frame'])*percent)):
        del_xy = (pivot - pivot.shift(tau+1, axis=0))**2
        en_msd.append((del_xy['x pos'] + del_xy['y pos']).stack().mean())
    return en_msd

#%pdb on
temp_track = tracked_data[0][0]
temp_track = temp_track.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
temp_track['track id'] = temp_track['track id'].astype('int64')
temp_track_sorted = temp_track.sort_values(['frame', 'track id'])
#temp_track.groupby('track id').apply(lambda x: x.name)

Calculated the ensemble averaged MSD on the first 3rd of the data.

'''Calculate the ensemble MSD for each window'''
en_msd = []
for tracks, window in tracked_data:
    temp_track = tracks.drop('track id', axis=1)
    temp_track = temp_track.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
    en_msd.append([calc_en_msd(temp_track),window])
    print 'done with window '+str(window)

done with window 3
done with window 4
done with window 5
done with window 6
done with window 7

'''Checkpoint after calculating ensemble MSD'''
import cPickle
pik = open('Ana_16042401_ensemble_msd.pkl', 'w')
cPickle.dump(en_msd, pik)
pik.close()

import cPickle
pik = open('Ana_16042401_ensemble_msd.pkl', 'r')
en_msd = cPickle.load(pik)
pik.close()

for num, entry in enumerate(tracked_data):
    df, window = entry
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d(df['x pos']%1, df['y pos']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    
    
    len_msd = len(en_msd[num][0])+1
    ax2 = plt.subplot(122)
    plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd), en_msd[num][0])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    
    plt.show()

SPIFF Correction

Do the SPIFF correction on the Mosaic Data

def spiff_correction(df, x_key='x', y_key='y'):
    df_copy = df.copy()
    spiff = df_copy[[x_key, y_key]] % 1
    spiff['spiff_x_corr'] = spiff[x_key].rank()/len(spiff)
    spiff['spiff_y_corr'] = spiff[y_key].rank()/len(spiff)
    df_copy[x_key] = np.floor(df_copy[x_key]) + spiff.spiff_x_corr
    df_copy[y_key] = np.floor(df_copy[y_key]) + spiff.spiff_y_corr
    return df_copy

corr_tracked_data = []
for tracks, window in tracked_data:
    tracks_temp = spiff_correction(tracks, x_key='x pos', y_key='y pos')
    corr_tracked_data.append([tracks_temp, window])

'''Calculate the ensemble MSD for each window (for only first third of tau)'''
corr_en_msd = []
for tracks, window in corr_tracked_data:
    temp_track = tracks.drop('track id', axis=1)
    temp_track = temp_track.rename(columns={'x': 'x pos', 'y': 'y pos', 'particle': 'track id' })
    corr_en_msd.append([calc_en_msd(temp_track),window])
    print 'done with window '+str(window)

done with window 3
done with window 4
done with window 5
done with window 6
done with window 7

MSD of Corrected SPIFF

for num, entry in enumerate(corr_tracked_data):
    df, window = entry
    plt.figure(figsize=[10,4])
    plt.subplot(121)
    plt.hist2d(df['x pos']%1, df['y pos']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(corr_en_msd[num][0])+1
    ax2 = plt.subplot(122)
    plt.plot(range(1,len_msd), corr_en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], corr_en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    print 'Average particles per frame = '+str(df.groupby('frame').apply(len).mean())
    
    plt.show()

Average particles per frame = 4338.1025641

Average particles per frame = 4357.97333333

Average particles per frame = 4313.88512821

Average particles per frame = 4367.75179487

Average particles per frame = 4301.75179487

Summary of SPIFF Correction and MSD

for num, entry in enumerate(corr_tracked_data):
    df, window = tracked_data[num]
    plt.figure(figsize=[8,6])
    plt.subplot(221)
    plt.hist2d(df['x pos']%1, df['y pos']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(en_msd[num][0])+1
    ax2 = plt.subplot(222)
    plt.plot(range(1,len_msd), en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax2.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax2.transAxes)
    plt.title('MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    df, window = entry
    plt.subplot(223)
    plt.hist2d(df['x pos']%1, df['y pos']%1, bins=50)
    plt.axis('equal')
    cbar = plt.colorbar()
    cbar.set_label('counts')
    plt.minorticks_on()
    plt.title('Corrected SPIFF for window='+str(window))
    plt.ylabel('Meta Pixel Y')
    plt.xlabel('Meta Pixel X')
    
    len_msd = len(corr_en_msd[num][0])+1
    ax4 = plt.subplot(224)
    plt.plot(range(1,len_msd), corr_en_msd[num][0], '-ob')
    pwr_law = lambda x, param: param*x
    fit_params = scipy.optimize.curve_fit(pwr_law, np.arange(1,len_msd)[-200:], corr_en_msd[num][0][-200:])
    plt.plot(np.arange(1,len_msd), pwr_law(np.arange(1,len_msd), fit_params[0][0]), 'g')
    plt.loglog()
    ax4.text(0.5, 0.2, str(fit_params[0][0])+'t', transform=ax4.transAxes)
    plt.title('Corrected MSD and Power Law for window='+str(window))
    plt.ylabel('MSD (Pixels$^2$)')
    plt.xlabel('Lag Time (Frames)')
    plt.tight_layout()
    
    print 'Average particles per frame = '+str(df.groupby('frame').apply(len).mean())
    
    plt.show()

Average particles per frame = 4338.1025641

Average particles per frame = 4357.97333333

Average particles per frame = 4313.88512821

Average particles per frame = 4367.75179487

Average particles per frame = 4301.75179487